library(tidyverse)
library(lfe)

rm(list= ls())

treat_date <- as.Date('2016-08-06')

## Load BA data

ba <- readRDS('data/BA_employment_muni.rds') %>% 
  filter(date > as.Date('2010-01-31')) %>%
  dplyr::select(date, treated, matches('emp_ref'),
                emp_nat_pc,
                emp_ref_pc,
                state, county, AA.district.id,
                date.id, ags) %>% 
  mutate(post = ifelse(date >= treat_date, 1, 0)) %>% 
  mutate(period = as.numeric(factor(date))) %>% 
  mutate(treated_post = treated * post) %>%
  filter(state %in% c('Bayern', 'Nordrhein-Westfalen'))

## Nr of pre and post treatment periods

id_treat <- ba %>% 
  filter(post == 1) %>% 
  filter(period == min(period)) %>% 
  pull(period) %>% unique()

ba <- ba %>% 
  filter(period >= (id_treat-11))

## Create Period*Treated Indicator Variables 

ba <- fastDummies::dummy_cols(ba, select_columns = 'period') %>% 
  mutate_at(vars(contains('period_')), ~ .*treated) 

## Generate Model Formulas 

outcomevars <- names(ba)[str_detect(names(ba), 'emp')]
periodvars <- names(ba)[str_detect(names(ba), 'period_')]

## Exclude last pre-treatment period 

periodvars <- periodvars[-which(periodvars == 'period_26')]

results <- lapply(outcomevars, function(dv){
  
  model_formula <- as.formula(paste0(dv, ' ~', paste0(periodvars, collapse = '+'), 
                                     '| period + ags |0 | AA.district.id'))
  
  model_est <- felm(model_formula, data = ba)
  n <- model_est$residuals %>% na.omit() %>% length()
  
  model_est <- model_est %>%
    broom::tidy(conf.int = T) %>% 
    filter(str_detect(term, 'period')) %>% 
    mutate(outcome = paste(dv)) %>% 
    mutate(n = n)
  
}) %>% reduce(rbind) 

## Rename the period variables

results <- results %>% 
  mutate(period = parse_number(term) - (id_treat - 1))

## Subset to per-capita measures 

results <- results %>% 
  filter(str_detect(outcome, 'pc')) %>%
  mutate(outcome_label = recode(outcome,
                                `emp_ref_pc` = 'Refugee Employment\n (per 1,000 capita)',
                                `emp_nat_pc` = 'Overall Employment\n(excluding refugees)')) %>%
  filter(period >= -10)

## Plot this

p1 <- ggplot(results, aes(x = period, y = estimate, 
                          ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0, linetype = 'dotted', color = 'grey40') +
  geom_vline(xintercept = 0, linetype = 'dotted', color = 'grey40') +
  geom_errorbar(width = 0) +
  geom_point(fill = 'white', shape = 21) +
  theme_bw() + 
  ylab('Coefficient Estimate') + 
  xlab('Quarter relative to treatment ') + 
  facet_wrap(~ outcome_label, scales = 'free') + 
  scale_color_grey(name = '') +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
  scale_x_continuous(breaks = seq(min(results$period), max(results$period), 1))

p1

#### Table A.2 - details on the analysis ####

l_out <- results %>% 
  dplyr::select(period, outcome_label, estimate, 
                std.error, p.value, n) %>% 
  mutate(across(3:6, ~round(., 3))) %>% 
  mutate(n = base::format(n,big.mark = ',')) %>% 
  mutate(outcome_label = str_replace_all(outcome_label, '\n', ' '))

## Kable 

library(kableExtra)

l_out %>%
  kableExtra::kable(., format = 'latex', 
                    booktabs = T,
                    linesep = "",
                    caption = 'Details on employment results',
                    col.names = c('Period', 'Outcome',
                                  'Estimate', 'SE', 'p-value', 'N')) %>%
  kable_styling(latex_options = c("hold_position", 'scale_down'), 
                position = "center") 
